function [paramsOpt,theta2Adj,input] = runInduceStatVAR_threshold_macro_spec(theta2,...
    VarU_x2,Cov10U_x2,Cov01U_x2,x2,macros,econAct,thresholdLevel,switch_h0,switch_hx)
dispOn = 1;
factors     = [macros;x2];
[nx,T]      = size(factors);
nm          = size(macros,1);
nx2         = size(x2,1);
VarU        = zeros(nx,nx,T);
Cov10U      = zeros(nx,nx,T);
Cov01U      = zeros(nx,nx,T);
for t=1:T
    VarU(nm+1:nm+nx2,nm+1:nm+nx2,t)   = VarU_x2(:,:,t);
    Cov10U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov10U_x2(:,:,t);
    Cov01U(nm+1:nm+nx2,nm+1:nm+nx2,t) = Cov01U_x2(:,:,t);
end

if any(switch_h0 == 1)
    h0RegimeOn = 1;
else
    h0RegimeOn = 0;
end
if any(any(switch_hx == 1))
    hxRegimeOn = 1;
else
    hxRegimeOn = 0;
end
[~,~,hx_1,hx_2] = unfoldTheta2_threshold_macro(theta2,size(factors,1));
if hxRegimeOn == 1
    % Regime 1
    eigValues_1         = eig(hx_1);
    modulus_1           = sqrt(real(eigValues_1).^2 + imag(eigValues_1).^2);
    if any(modulus_1>=1)
        params.delta1 = 1/max(modulus_1)*0.99;
    else
        calibrateParams.delta1 = 1;
    end
    % Regime 2
    eigValues_2         = eig(hx_2);
    modulus_2           = sqrt(real(eigValues_2).^2 + imag(eigValues_2).^2);
    if any(modulus_2>=1)
        params.delta2 = 1/max(modulus_2)*0.99;
    else
        calibrateParams.delta2 = 1;
    end
else
    % We have no regime-dependent loadings, so we just check stability
    % using hx_1
    eigValues_1         = eig(hx_1);
    modulus_1           = sqrt(real(eigValues_1).^2 + imag(eigValues_1).^2);
    modulus_2           = 0; 
    if any(modulus_1>=1)
        params.delta1 = 1/max(modulus_1)*0.99;
    else
        calibrateParams.delta1 = 1;
    end
    calibrateParams.delta2 = 1;
end
if ~exist('calibrateParams','var')
    calibrateParams = [];
end

if any(modulus_1>=1) || any(modulus_2>=1) 
    if dispOn == 1
        disp('OPTIMIZATION TO INDUCE STATIONARITY OF THRESHOLD MODEL')
    end
    % Moment-based scaling to induce stationarity
    meanFactors           = mean(factors,2);
    T                     = size(factors,2);
    input.varXData        = (factors-repmat(meanFactors,1,T))*(factors-repmat(meanFactors,1,T))'/(T-1)-mean(VarU,3);
    input.VarU            = VarU;
    input.Cov10U          = Cov10U;
    input.Cov01U          = Cov01U;
    input.factors         = factors;
    input.theta2          = theta2;
    input.selectParams    = fieldnames(params);
    input.calibrateParams = calibrateParams;
    input.econAct         = econAct;
    input.thresholdLevel  = thresholdLevel;
    input.switch_h0       = switch_h0;
    input.switch_hx       = switch_hx;
    input.h0RegimeOn      = h0RegimeOn;
    input.hxRegimeOn      = hxRegimeOn;
    input.numSim          = 1D5;
    input.burnin          = 1D4;
    maxIter               = 10000;
    tolFun                = 1e-6;
    options       = optimset('Display','iter','MaxIter',maxIter,'MaxFunEvals',10000,'TolFun',tolFun,'TolX',1e-6);
    paramsValues  = struct2array(params)';
    
    paramsValuesOpt = fminsearch(@induceStatVAR_threshold_macro_spec,paramsValues,options,input);
    [Qopt,theta2Adj,paramsOpt,jointDownScalingOn] = induceStatVAR_threshold_macro_spec(paramsValuesOpt,input);
    if jointDownScalingOn == 1
        disp(['Joint downscaling of hx_1 and hx_2 by ', num2str(min(paramsValuesOpt))])
        paramsOpt.delta1 = min(paramsOpt.delta1,paramsOpt.delta2);
        paramsOpt.delta2 = min(paramsOpt.delta1,paramsOpt.delta2);
    else
        for i=1:size(paramsValuesOpt,1)
            name = input.selectParams{i};
            disp([name ' = ', num2str(paramsValuesOpt(i,1)), ' to induce a stationary threshold model. Qopt = ', num2str(Qopt)]);
        end
    end
else
    paramsOpt.delta1 = 1;
    paramsOpt.delta2 = 1;
    theta2Adj        = theta2;
end

end